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This study explores a method to characterize temporal structure of intermittent phase locking 
in oscillatory systems. When an oscillatory system is in a weakly synchronized regime away from 
a synchronization threshold, it spends most of the time in parts of its phase space away from 
synchronization state. Therefore characteristics of dynamics near this state (such as its stability 
properties/Lyapunov exponents or distributions of the durations of synchronized episodes) do not 
describe system's dynamics for most of the time. We consider an approach to characterize the system 
dynamics in this case, by exploring the relationship between the phases on each cycle of oscillations. 
If some overall level of phase locking is present, one can quantify when and for how long phase 
locking is lost, and how the system returns back to the phase-locked state. We consider several 
examples to illustrate this approach: coupled skewed tent maps, which stability can be evaluated 
analytically, coupled Rossler and Lorenz oscillators, undergoing through different intermittencies 
on the way to phase synchronization, and a more complex example of coupled neurons. We show 
that the obtained measures can describe the differences in the dynamics and temporal structure of 
synchronization/desynchronization events for the systems with similar overall level of phase locking 
and similar stability of synchronized state. 
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PACS numbers: 05.45.Xt, 05.45.Tp, 05.45.-a 



I. INTRODUCTION 



Synchronization of oscillations is a widespread phe- 
nomenon in physical sciences and beyond [3,13], includ- 
ing physiology J|| and neuroscience 0, @ > and has being 
studied using approaches and methods of physics and 
nonlinear dynamics even in that applied context [(| Q . In 
particular, phase synchronization, dynamical phenomena 
where the phases of oscillations approach each other with 
time while both phases may remain chaotic and their am- 
plitudes remain uncorrelated, appears to be quite com- 
mon in nature and is well-studied [3]. It is known that 
for some parameter values (usually, moderate coupling 
strength), phase synchronization can occur in an inter- 
mittent fashion |8-ll0|. It will present itself as an in- 
termittent phase locking. The phase difference of two 
oscillators is close to constant during finite time inter- 
vals. These intervals of temporal synchronization epochs 
are random and interspersed by desynchronization events 
which arc characterized by varying phase differences. Dif- 
ferent types of intermittent behaviors may take place in 
different systems fioj and transition through intermit- 
tency may depend on the geometric structure of chaotic 
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attractors [13] ■ Three types of intermittencies have been 
observed near the onset of phase synchronization: the 
type-I intermittency [|[ , the eyelet intermittency [1, [13] , 
and the ring intermittency @ . 

Different approaches to synchronization analysis exist. 
One may approach the problem from the time-series anal- 
ysis perspective, characterizing how much time-series are 
correlated in some particular sense. This view essentially 
tells how close we are to the synchronized state. Another 
approach is to study the stability of the synchronized 
state, via, for example, Lyapunov exponents, so that 
small positive Lyapunov exponents will describe the weak 
instability and, potentially, intermittency of synchroniza- 
tion. Lyapunov exponents, distribution of durations of 
synchronized episodes etc. are ways to explore the syn- 
chronized state/synchronization manifold. Certain uni- 
versality of behavior of a weekly unstable synchronization 
is expected and is captured by these approaches. Studies 
of Lyapunov exponents spectra [l3j may go further and 
explore the synchronous/desynchronous states dynamics 
in the vicinity of the saddle- node bifurcation point. How- 
ever, these approaches become less powerful when the 
system is far away from bifurcation point and do not 
describe the desynchronized episodes of intermittent or 
otherwise variable and weak synchronization in a wide 
range of parameter space. The desynchronized dynamics 
("reinjection mechanisms" in intermittency terminology) 
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may be not universal and depend on the peculiarities of 
the system under study. 

If the system is in imperfectly synchronized state, but 
is still close to synchronization threshold (the case, which 
is probably expected in many physics or engineering ap- 
plications), it spends most of the time in an almost 
synchronized state (near synchronization manifold) and 
the focus on this state is natural. However, the syn- 
chrony can be weak. In particular, in living systems 
the synchronization may be very imperfect and weak, 
and too strong synchrony (high degree of temporal co- 
ordination) may be responsible for pathological states, 
such as schizophrenia or Parkinson's disease [3 EH- In 
terms of the phase space, the system spends relatively 
small fraction of time near synchronization manifold in 
the phase space. Thus even if Lyapunov exponents or 
distributions of the durations of synchronized episodes 
may be informative about this part of dynamics, they 
do not tell much about dynamics overall, because it is 
mostly desynchronized. This calls for the study of the 
highly- variable synchrony and for the study of the struc- 
ture of the phase space away from the relatively strongly 
unstable synchronization manifold. One may not expect 
much of universality here, but the characterization of this 
dynamics should be possible. 

One approach is to characterize the degree of syn- 
chronization for a series of sliding time-windows and ob- 
tain statistical estimates of significance (as in. for exam- 
ple, [lH ). This may reveal important details of synchro- 
nized dynamics, not accessible otherwise (as in [It} for 
the case of tremor oscillations). However, synchroniza- 
tion is not instantaneous phenomena, so this approach is 
not designed to detect changes in the dynamics on the 
very short time-scales. Yet, changes of this kind (very in- 
termittent, variable synchronization) have been observed 
experimentally and have been conjectured to have func- 
tional significance [HI, ■ 

Another example, where the fine temporal details of 
synchrony may be relevant may come from ecological dy- 
namics. Synchronization in ecological dynamics is hardly 
perfect. It was conjectured that prolong zero-lag synchro- 
nization elevates risks of extinction, so that there is an 
interest in these cases (reviewed in, e.g. 20]). One may 
want to distinguish between the moderate level of syn- 
chrony with few prolong synchronized episodes and many 
short synchronized episodes. 

Here we expand on the approach of [IH and study 
a method to analyze intermittent and weak phase lock- 
ing from the nonlinear dynamics perspective. This ap- 
proach is based on the constructions of the first-return 
maps for the difference of phases, partition of the result- 
ing phase space into synchronized and non-synchronized 
regions and the analysis of transitions between them. We 
show that for the case when some overall level of phase 
locking is present, one can characterize how this phase 
locking (properly defined) is gained, maintained or lost 
for each cycle of oscillations and what it means in terms 
of the organization for the phase space of coupled oscil- 



lators system. We do so with examples drawn from sev- 
eral typical model systems, ranging from simple coupled 
tent maps to more realistic conductance-based models of 
coupled cells. We also show how one can characterize the 
differences in the dynamics of coupled systems, which un- 
dergo through the same kind of intermittency and posses 
the same stability properties of synchronization manifold. 

The structure of this paper is the following. Section [TT1 
presents the set up of the method. Section Hill considers 
the method as applied to the dynamics of simple coupled 
maps. Section IIVI considers different types of intermit- 
tencies for phase synchronization. Section [V] illustrates 
the ideas of the paper for a more complex dynamical sys- 
tem. Finally, in Section IVII we discuss the main ideas 
behind the method, summarize its properties, and dis- 
cuss its applicability and possible further developments. 

II. METHODS TO STUDY THE TEMPORAL 
VARIABILITY OF THE SYNCHRONIZATION 

Suppose that two oscillatory signals allow for a reason- 
ably good reconstruction of their phases (f>i (t) and cj>2 (t) . 
In a discrete-time version (and thus applicable to exper- 
imentally obtained data) one can consider a standard 
index to characterize the strength of the phase locking 
between these two signals: 

7=II^E^ (tj) H 2 > W 
i=i 

where 0(tj) = (fti{tj) — 4>2{tj) is the phase difference. 
This synchronization index varies from to 1, that is 
from complete lack of phase locking to perfect phase lock- 
ing (and does not tell anything about amplitudes of the 
signals). One may also compute Lyapunov exponents, 
which characterize stability of the synchronization man- 
ifold. In general, negative values of all transversal Lya- 
punov exponents are not necessarily sufficient to guaran- 
tee perfect synchronization nevertheless their values 
describe the degree of stability/instability of the synchro- 
nized state and the strength of synchronization. 

As we explained in the Introduction, these measures 
are not designed to describe the fine temporal structure 
of the dynamics, which motivates the following approach 
(its major steps were sketched by us in [l8[ and its expla- 
nation, validation and properties are presented in subse- 
quent sections of this paper). First, one needs to confirm 
that there exist some degree of phase locking between the 
phases of two variables x and y. This can be detected, for 
example, by computing the Eq. (TTJ) and comparing the re- 
sulting value of 7 against some kind of significance value, 
computed analytically or via surrogates [l6[ . It should be 
emphasized here that our method allows for an analysis 
of the temporal development of phase difference if some 
level of synchrony (some preferred phase-locking angle) 
is present. It makes no sense if there were no synchrony 
between the signals at all. 
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Then set up a check point for the phase of y: when- 
ever the phase of y crosses this check point from below, 
record the phase value of x, generating a set of consec- 
utive phase values {4> x ,i}, i = 1, . . . , iV, where N is the 
number of such level crossings in an episode under anal- 
ysis. Then plot <j) x .i+i vs. (f> x .i for i = I, . . . , N — 1. The 
properties of these plots will characterize the dynamics 
of the synchronization. A tendency for predominantly 
synchronous dynamics will appear as a cluster of points, 
with the center at the diagonal 4> x ,i+i = <f>x,i- 

After determining the center of the cluster (that is, 
determining the preferred phase difference between two 
signals, at which they have the tendency to be locked) 
all values of the phases are shifted to a position with the 
center of the first region (quadrant) (see Fig. [TJ. This 
step is not necessary and is done here for the uniformity of 
the analysis. Then we consider how the system leaves the 
synchronization cluster and its vicinity and how it returns 
back to the synchronization by quantifying transitions 
between different regions of the (4> x ,i, 4>x,i+i) space. 
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FIG. 1. Diagram of the (cj> x ,i, <j>x,i+i) first-return map. The 
arrows indicate all possible transitions from one region to an- 
other and the expressions next to the arrows indicate the rates 
for these transitions. 

This phase space is then partitioned into four regions 
as illustrated in Fig. Q] (let us remind, that the space of 
) is in fact torus). While the system is in 
the first region, it is considered to be in a synchronized 
state. Dynamics outside of the first region will be called 
a desynchronized state. Thus the synchronized state here 
is the one, where the deviation from the preferred phase 
angle is less than 7r/2. This is a compromise value: not 
too small to allow for some moderate fluctuations of the 
phases; not too large to allow for the phases to be suf- 
ficiently related and to be involved in some function. 
It also allows for symmetric partitioning of the phase 
space and definition of a few easy-to-compute transition 
rates (see below). Some situations may require a dif- 
ferent definition of the synchronized state and the rest 



of the method will need to be modified. Nevertheless, 
even though this partition is quite coarse, it allows for 
an inspection of the fine temporal details, as we will 
show below. Four resulting regions in the phase space 
are numbered in a clockwise manner (Fig. [1} , since this 
is the primary direction of the dynamics. An example of 
the first-return map for two coupled Lorenz oscillators is 
presented in Fig. [5] (detailed study of this system with the 
method described here is in the Section UVl sec Fig. [7]). 

Transition rates {i = 1, 2, 3, 4) for the transitions be- 
tween four regions of the map are defined as the number 
of points in a region, from which the system leaves the re- 
gion, divided by the total number of points in that region. 
For example, n is the ratio of the number of trajecto- 
ries escaping the first region for the second region to the 
number of all points in the first region. While the time- 
averaged measures of synchrony characterize whether the 
synchronization is strong or weak overall, utilization of 
these rates lets us explore the dynamics of the synchro- 
nization in time. 
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FIG. 2. An example of the first-return map for two coupled 
Lorenz oscillators. All four first-return plots have the same 
data points (gray circles), but each subplot (a-d) presents 
the evolution of points from one region (I, II, III, and IV, 
respectively). If a point evolves from one region to another 
region, then we represent it as o — *. If a point evolves within 
the same region, then we represent it as o — o. Thus each plot 
shows the transitions from a corresponding part of the phase 
space. We can compute the transition rates ri,2,3,4 from the 
panels (b), (d), (c), and (a), respectively. 

Further exploring the properties of the dynamics in 
the space (4>x,it <Am+i) ; one ma y compute the relative 
frequencies of desynchronization events of different dura- 
tions. In the considered first-return map approach, the 
duration of a desynchronization event is the number of 
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steps that system spends away from the first region mi- 
nus one (because, the point on the map has two coordi- 
nates, one of which is a "future" phase). The shortest 
duration of the desynchronization event corresponds to 
the shortest path 2-4-1 (note that the desynchronization 
will always start at the second region by the virtue of 
our description of the dynamics). This will correspond 
to the desynchronization length of one cycle (in other 
words, in two cycles the phases are back in a locked 
state). If we want to consider the chance of the dura- 
tion of two cycles, we should consider the path 2-3-4-1. 
Longer desynchronization events will have many differ- 
ent paths corresponding to them. We will systematically 
compute the relative frequencies of the desynchroniza- 
tion events of different durations in the Section IIVI for 
different systems. 

The transition rates r2,3,i are related to the durations 
of desynchronizations. Higher values of ^2,3,4 promote 
short desynchronization episodes, while their low val- 
ues promote long desynchronization events (note that 
for a well-synchronized dynamics, there may be very few 
points in the third quadrant, and numerical estimation 
of r% may be unreliable). If transitions are independent 
then the relative frequencies of the durations of desyn- 
chronization events can be estimated by multiplying the 
corresponding rates of the transitions for different paths. 
For example, to estimate the probability of a desynchro- 
nization event of the shortest duration we should consider 
the shortest path 2-4-1 and the corresponding probabil- 
ity is r2r4. This will correspond to the desynchronization 
length of one cycle. Longer desynchronizations will re- 
quire summation of the products of rates corresponding 
to longer paths of equal length. In general, this inde- 
pendence (or near independence) may or may not be in 
place, but we consider several examples, where the dif- 
ference is not large and the rates ^,3,4 really describe the 
desynchronization durations. 

The rate r\ is related to the averaged duration of the 
laminar phase and thus characterizes the property of syn- 
chronized state. The average duration of the synchro- 
nized phase < £ > (for an intermittency near onset of 
phase synchronization, this quantity usually scales in a 
universal manner for specific intermittency types) is in- 
versely proportional to r\ : 

<*>~l/n (2) 

and < £ >= l/r±, if < £ > is measured in the number of 
iterations of the map (j> x i, which is essentially the number 
of cycles of oscillations. The closer r\ is to 1, the more 
frequently the synchronized dynamics is interrupted. 

Thus the rates permit evaluation of whether the weak- 
ness of synchrony is due to a few relatively long desyn- 
chronization events or to large number of short desyn- 
chronization events. 



III. COUPLED MAPS EXAMPLE 

In this section we consider behavior of rates for dif- 
ferent parameter values (and thus for different transver- 
sal Lyapunov exponents of synchronized state and the 
synchronization index 7) in a simple system. We con- 
sider linearly coupled skew tent maps. This relatively 
simple map allows for exact computation of Lyapunov 
exponents. The properties of the phase space depend 
on the parameters in a simple way. Thus we are able 
to see the relationship between system's parameters, dy- 
namics, Lyapunov exponents and the rates r*j. This sys- 
tem is not necessarily an appropriate choice for the study 
of phase synchronization as a dynamical phenomenon, 
when phases are correlated and amplitude are not plj . 
However we are not as much concerned with the phase 
synchronization per se in this sense, as we are concerned 
with a simple treatable example of synchronous dynamics 
here. Therefore, coupled skew tent maps prove a possi- 
bility to illustrate our ideas in a simple system, where 
analytical treatment is partially possible. 

Consider a skew tent map 

'<-*>- ft- <3) 

where < a < 1. We will use two such maps with vari- 
ables x and y. Let us consider the following coupling: 

x(t + 1) = (1 - e)f(a, x(t)) + ef{a, y(t)), (4) 
y(t + 1) = ef(a, x(t)) + (1 - e)f(a, y(t)), 

where e is the coupling strength. 

To characterize synchrony between x and y, we con- 
sider a new variable 

e(t)=y(t)-x(t), (5) 

which will be an analog of the phase difference described 
above. In these coupled maps, synchronous state x = y 
is stable for e > e c where the critical value e c < 1/2. 
Note that 9 in Eq. varies between —1 and 1 while the 
phase variable <pi m Fig. [1] varies between — 7r and 7r. In 
this section, we rescale the phase diagram accordingly to 
define the transition rates r-s. 

Two Lyapunov exponents can be computed analyti- 
cally (see e.g. [l|): 

A(a) = — olno — (1 — a) ln(l — a), (6) 
Aj_(a, e) = —a In a — (1 — a) ln(l — a) + In |1 — 2e|. 

The first Lyapunov exponent A(a) corresponds to the mo- 
tion governed by the one-dimensional dynamical system 
x(t + 1) = f(a,x(t)) whose subspace is {(x,y) | x = 
y}. The second, transversal Lyapunov exponent X±(a,e) 
characterizes the dynamics transverse to that invariant 
subspace. Note that A(a) is always positive for < a < 1 
and has a maximum In 2 at a = 1/2. Since Aj_(a,e) = 
A (a) + In 1 1 — 2e\, the transversal Lyapunov exponent 
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X±(a,e) can change from positive to negative depend- 
ing on a and s. Both Lyapunov exponents are symmetric 
around a = 1/2. That is, A(a) = A(l — a) and Aj_(a, s) = 
Aj_(l — a,e). The synchronization stability threshold 
Aj_(a,e c ) = is reached at e c = | — ^a a (l - a)^ 1_a ^. 

It is interesting to note that the coupled maps (Eqs. ([3]) 
and ([3])) have exact solutions when the system is at a 
bifurcation point and is not robust. For example, for with 
a = 1/2 and e = 1/4 (and thus A = In 2 and A± = 0) 
initial conditions (x,y) = (0.46,0.5) lead to a period-20 
trajectory and initial conditions (x,y) = (0.044,0.052) 
lead to a period-100 trajectory. These trajectories are 
not hard to compute, and then one can formally compute 
Ti (ri = 2/5, r 2 = 3/5, r 3 = 2/5, and r± = 2/5 for the 
former trajectory and r\ = 13/25, r 2 = 12/25, r 3 = 
13/25, and r± = 13/25 for the latter). However, the 
system docs not exhibit synchronous behavior, there is 
no single preferred value for the difference of variables 9, 
and computation of rates is not appropriate. 

We numerically iterate the coupled maps (Eqs. J3]) 
and ((¥])) and compute the rates for different a and 
e. We vary a from 0.125 to 0.35 with step size 0.025 and 
e = (1 — ka a (\ — a) 1 ~ a )/2 where k varies from 1.05 to 
1.45 with step size 0.05 (with this set-up A^ = lnfc). 
Fig.[3]shows the relationship between Aj_ and the rates r;. 
Each dot represents a different pair of (a,e). Since Aj_ 
characterizes the evolution of the trajectory transverse 
to the invariant subspace {(x,y) \ x = y}, the smaller 
Aj^ corresponds to less unstable manifold and more syn- 
chronized dynamics, which results in the smaller r\. As 
Aj_ decreases, the points (x, y) coalesce to the diagonal 
line x = y in the (x,y) plane. This shortens the dura- 
tion of desynchronization events and increases the rein- 
jection probability to the synchronous state. Thus, r\ 
tends to decrease while r 2 ,4 tend to increase as A^ de- 
creases. However, one can clearly sec that the same value 
of Aj^ may correspond to drastically different rates and 
thus to different timing of synchronized/desynchronized 
dynamics. 

If there is no point in the given region, then the tran- 
sition rate at that region is undefined. As X± decreases, 
the system spends less and less time in desynchronous 
regions. This makes the number of points in the third 
region extremely small (or sometimes zero). Thus nu- 
merically computed values of r 3 are less reliable in this 
case (or are undefined, as is the case of r 3 for small Aj_ 
in Fig. He)). 

We now consider some specific parameter values to 
show that the transition rates can effectively discrimi- 
nate two different systems with identical Ax and simi- 
lar 7. Fig.[4ja) shows coupled skew tent maps for a = 0.3 
and a = 0.7. Their two Lyapunov exponents are exactly 
the same because of their symmetry around a = 1/2, 
but the maps are clearly different and their dynamics are 
different. Fig. SJb) illustrates how coinciding transver- 
sal Lyapunov exponents Ax (stability property) and syn- 
chronization indices 7 (average strength of phase locking) 
for both maps depend on e. For these values of a the 
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FIG. 3. We varied parameters a and e to compute the 
transversal Lyapunov exponent Ax and rates ri. The rate 
ri (r-2) and Ax have a positive (negative) relation. Rates 
and ta present more spread and more complex relationship 
with Ax- 



threshold value of coupling is e c « 0.2286. 
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FIG. 4. (Color online) (a) Two skew tent maps for a = 0.3 
(thick solid line) and a — 0.7 (thin solid line), respectively, 
(b) Transversal Lyapunov exponent Ax (dash-dotted line) and 
two synchronization indices 7 (thick solid line for a — 0.3 and 
thin solid line for a = 0.7). Two 7 values are almost the same. 
Transversal Lyapunov exponent Ax decreases as e increases 
and crosses 0, when 7 approaches 1. 

However, Fig. [S] shows that the rates r, for a = 0.3 
and a = 0.7 are different and sometimes substantially 
different for a range of e € (0, 0.15). The values of the rate 
r\ are essentially the same, as one may expect, because 
T\ is related to the stability of synchronized state, which 
is the same for both cases. The other rates present a 
different picture. The system with a = 0.7 (thin solid 
line in Fig. [S]) has higher r2. 3j 4 for almost all range of e 
than those of a = 0.3 (thick solid line), which suggests the 
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stronger tendency for shorter desynchronization events. 

As the coupling e increases to synchronization thresh- 
old value, the rates ^2,3,4 are undefined while r± decreases 
to zero. This is expected because the system spends most 
of (if not all) the time in the vicinity of the synchronized 
manifold. Thus, the dynamics is mostly synchronized 
and the rates are not much relevant. The rate r± sa at 
e sw 0.15 < e c « 0.2286. This probably happens because 
our partition of the phase space of the first-return map 
for the phase difference implies that r% = and dynam- 
ics is synchronous when the deviation from the complete 
synchrony 8 = is less than 1/2. But 9(t) can still fluc- 
tuate in some smaller limits, which is what happens for 
e in the interval of about (0.15,0.2286). 
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FIG. 5. (Color online) Transition rates for the linearly coupled 
skew tent maps for a = 0.3 (thick solid line) and a — 0.7 (thin 
solid line). We used 350000 iterations and removed the first 
50000 iterations. Differences in the rates ^"2,3,4 are clearly 
visible. 

This simple case illustrates that although the stabil- 
ity of the synchronization manifold is the same for two 
systems, their dynamics can be substantially different. 
Lyapunov exponents and synchronization index 7 natu- 
rally arc unable to capture the differences. However, the 
analysis using the transition rates can distinguish two 
different systems and provide the temporal details of dy- 
namical behaviors. 



IV. ANALYSIS OF THREE DIFFERENT 
INTERMITTENCY TYPES 



In this section, we will study these three types of in- 
termittent behaviors in coupled chaotic oscillators with 
techniques described in Section [TTJ First, we will look 
at the type-I and the eyelet intermittencies in bidirec- 
tionally coupled Rossler oscillators and bidirectionally 
coupled Lorenz oscillators. And then, we will look at 
the eyelet and the ring intermittencies in unidirection- 
ally coupled Rossler oscillators. 

Intcrmittcncy types are characterized by the proba- 
bility distribution of the laminar (synchronized in our 
context) phase length £ and the scaling behavior of the 
mean laminar phase length < £ > . As far as the type- 
I and the eyelet intermittencies are concerned, coupled 
Rossler oscillators and coupled Lorenz oscillators with 
small frequency detuning show the same scaling prop- 
erties of synchronous episodes near and away from the 
point of bifurcation 0, |T3 - However, we will show that 
there are important differences in timing of dynamical 
behavior. 

We consider two coupled chaotic Rossler oscillators: 

^i, 2 = -^i, 22/i, 2 - 2=1,2 + £1, 2(2*2, 1 - 2:1,2), (7) 
2/1,2 = ^1,22:1,2 + 0.15yi, 2 , 
z' 12 = 0.2 + zi, 2 (a:i,2 - 10), 

where wj. 2 are the natural frequencies of each chaotic 
oscillator and £1,2 measure the strengths of the linear 
coupling. The phase difference is obtained by 

9 — (f>i — fa — arctan ( — ] — arctan ( — ] . (8) 
\xj \x 2 J 

We also consider two coupled chaotic Lorenz oscilla- 
tors: 

2'i,2 = 10-0(2/1,2 - 2:1,2) + £(2:2,1 - 2:1,2), (9) 
2/1,2 = (36.5 + 71,2)2:1,2 - 2/1,2 - £1,221,2, 

z[ 2 = -3.0zi, 2 + 2:1,22/1,2, 

where 71,2 are parameters for the frequency detuning of 
each chaotic oscillator and e measures the strength of 
the linear coupling. As in Q, by defining a new variable 
u = \/ x 2 + y 2 a phase can be defined on the [u, z) plane. 
The phase difference is obtained by 

e = 4>i - 02 (10) 

= arctan ( — ] — arctan ( — — ] , 

\U1-U-1J \U2-U2J 

where (u\, z±), (112, £2) are unstable fixed points in the 
middle of the trajectories rotating in the (it, z) plane. 



Intermittently synchronized dynamics can be observed 
near boundary of synchronization. So far, three types of 
intermittencies near the onset of phase synchronization 
have been observed: the type-I intermittency @ , the eye- 
let intermittency 0, [l2[ , and the ring intermittency . 



A. Type-I intermittency 

In coupled Rossler oscillators, let us assume that the 
difference between lo\ and W2 is small and the systems are 
bidirectionally coupled with £1.2 = £■ When uj\ = 1.015 
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and u>2 = 0.985, Lee et al. [8j found the first critical 
value e t = 0.0276 where 9 increases in a nearly periodic 
sequence of 2tt jumps. They showed that for e < e t , the 
mean laminar phase length < t > obeys the scaling rule 



(e t -e) 



-1/2 



(11) 



In coupled Lorcnz oscillators, let us again assume that the 
difference between 71 and 72 is small. When 71 = 1.5 and 
72 = —1.5, Lee. et al. Q also found the first critical value 
£t = 6.7. The phase desynchronization occurs with 2tt 
phase slip, but it shows ±27r irregular phase jumping be- 
haviors which arc different from coupled Rossler oscilla- 
tors. However, the probability distribution of the laminar 
phase length £ and the scaling rule of the mean laminar 
phase length < £ > are the same as coupled Rossler oscil- 
lators. The scaling rule together with the corresponding 
probability distributions of synchronous episodes suggest 
that these two systems are in the type-I intermittency. 
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FIG. 6. (Color online) Synchronized dynamics of bidirection- 
ally coupled Rossler oscillators in dependence on the cou- 
pling strength e. Two critical values for intermittencies are 
£t = 0.0276 and e c = 0.0286. (a) Rates Ti and synchroniza- 
tion index 7. (b) Relative frequencies of the desynchroniza- 
tion events of different durations (durations are measured in 
cycles of oscillations, i.e. in the number of iterations of the 
4> x ,i map, see Section [TT| ). Thick solid line represents the 
actual relative frequencies of the desynchronization events of 
different durations, while thin solid line represents the rela- 
tive frequencies, which would be there under assumption of 
independent transition rates. The longer cycles of desynchro- 
nization are dominant relative to short cycles. 

Figs. [5] and [7] show the rates ri and the relative fre- 
quencies of the desynchronization events of different du- 
rations for a range of coupling strength. By comparing 
Fig. [B] with Fig. one can tell that two coupled systems 
exhibit different temporal features of the dynamics. 

In Fig.^b), the probability of the long desynchroniza- 
tion events (longer than 5 cycles of oscillations) is in the 



range of about (0.8,0.9) while in Fig. [Jjb), this proba- 
bility is in the range of about (0,0.25). This is true, in 
particular, for similar values of 7. Since coupled Rossler 
oscillators have high probability of long cycles, two oscil- 
lators spend long time in the desynchronized state once 
the trajectory leaves synchronous state. However it does 
so very rarely. Coupled Lorenz oscillators (Fig.[7J present 
a different temporal dynamics. Short desynchronization 
cycles prevail and the comparable values of 7 are reached 
with a larger number of short desynchronized episodes. 
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FIG. 7. (Color online) Synchronized dynamics of bidirection- 
ally coupled Lorenz oscillators in dependence on the coupling 
strength e. Two critical values for intermittencies are et = 6.7 
and e c = 12. (a) Rates n and synchronization index 7. The 
behaviors of the rates and the durations of desynchronization 
events are different from coupled Rossler oscillators, (b) Rel- 
ative frequencies of the desynchronization events of different 
durations (durations are measured in cycles of oscillations, i.e. 
in the number of iterations of the c/> Xt i map, see Section HT|) . 
Thick solid line represents the actual relative frequencies of 
the desynchronization events of different durations, while thin 
solid line represents the relative frequencies, which would be 
there under assumption of independent transitions. Unlike 
Fig. [6l short desynchronization events dominate the dynam- 
ics. 



The relative frequencies of desynchronization events of 
different durations can be close to the relative frequen- 
cies estimates from the transition rates like Figs. [Tfb) 
and \5\b) or can be different from these estimates like 
Figs. [Bfb) andOJb). The transitions between regions are 
not necessarily independent (or near independent). How- 
ever, the estimated relative frequencies of desynchroniza- 
tion events of different durations still show some similar 
tendencies. 
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B. Eyelet inter mittency 

As the coupling strength increases further from the 
first critical value e t , exponentially rare 2tt phase jumps 
for coupled Rosslcr oscillators and irregular ±2tt phase 
jumps for coupled Lorenz oscillators can be observed be- 
fore the second critical value e c where the phase synchro- 
nization is established. Since the phase slips are rare, 
the scaling rule and probability distributions are differ- 
ent from those of the type-I intermittency. Lee et ah P 
found the second critical values e c = 0.0286 for coupled 
Rossler oscillators and e c = 1 2 for coupled Lorenz oscilla- 
tors. Note that other parameter values wi^ and 71.2 are 
the same as in the type-I intermittency. For £ t < e < s c , 
the mean laminar phase length < £ > obeys the scaling 
rule 



a r i r 2 r 3 r 4 



ln(<*>) ~ -(e c -e) 



1/2 



(12) 



The scaling rule and the corresponding probability distri- 
bution of the laminar phase length suggest both systems 
are in the eyelet intermittency. Since the eyelet intermit- 
tency features the longer laminar phase than that of the 
type-I intermittency, r\ is smaller than that of the type-I 
intermittency. 

Note that when e > £ c , the system is in the first region 
of the first-return map for phases, that is in the syn- 
chronous region, and the rates ^,3,4 are undefined while 
the rate r\ = 0. However the phase difference may fluc- 
tuate in some smaller limits and synchronization index 7 
may still be below 1. 

Figs. [DJa) and [3a) show that the transition rates Ti 
of both coupled Rossler oscillators and coupled Lorenz 
oscillators are characteristically distinct from one an- 
other near and away from the second critical value e c . 
In fact, the rate for coupled Rosslcr oscillators is sub- 
stantially smaller than that of coupled Lorenz oscillators 
for all range of e. This small value of r$ together with 
small value of S (0.25,0.4) for coupled Rossler oscil- 
lators implies that these two oscillators spend substan- 
tially longer time in the desynchronous state once the 
trajectory leaves synchronous state than coupled Lorenz 
oscillators do. The rates ^.3,4 for coupled Lorenz oscilla- 
tors show substantial variability while those for coupled 
Rossler oscillators do not. Therefore the temporal struc- 
ture of synchronized/desynchronized events is different 
in bidirectionally coupled Rossler oscillators and Lorenz 
oscillators even though the overall synchrony can be sim- 
ilar. 

Now let us also consider unidircctionally coupled 
Rossler oscillators for the eyelet intermittency. We use 
Eq. (JT)) with £\ — and £2 = e. When control parameters 
cji = 0.93 and U2 = 0.95, Hramov et al. [l2[ found two 
critical values £4 = 0.0345 and e c = 0.042. They showed 
that when e € (£t,£ c ), the intermittent behavior of uni- 
directionally coupled Rossler oscillators can be classified 
as the eyelet intermittency Note that this intermittent 
behavior can be also treated as the type-I intermittency 
with noise 12 1. 
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FIG. 8. (Color online) Eyelet intermittency in unidirec- 
tionally coupled Rossler oscillators with parameters ei = 0, 
£2 = £, and control parameters wi = 0.93, co2 = 0.95. Two 
critical points are et = 0.0345 and e c = 0.042. (a) Rates 
Ti and synchronization index 7. (b) Relative frequencies of 
desynchronization events of different durations. Thick solid 
line stands for actual duration, while thin solid line represents 
the estimates from the transition rates. The intermittent be- 
havior is different from the bidirectionally coupled Rossler 
oscillators. 



During the eyelet intermittency, rate r% in Fig. [Ha) 
experiences substantial variation while that in Fig. (Hlja) 
keeps almost the same level. Rate r± in Fig. [5Ja) also 
shows relatively large variation while that in Fig. (Hlja) 
does not. Note that rates r\ for both systems are almost 
the same. The large variations also can be observed in 
the durations of desynchronization events for several cy- 
cle lengths (1, 2, > 5) in Fig. [8jb) while almost no vari- 
ation in Fig. IHJb) is observed. Therefore, their temporal 
structures of synchronized/desynchronized events are dif- 
ferent even though the overall synchrony can be similar. 



C. Ring intermittency 

For the ring intermittency, we consider unidirectionally 
coupled Rossler oscillators. As in we use Eq. ([7]) with 
£1 = 0, £2 = £, and control parameters u>i = 1.0, LU2 = 
0.95. The frequency mismatch between two oscillators is 
larger than those of the type-I and the eyelet intcrmit- 
tencics discussed in the previous sections. To explore the 
synchronization behavior more clearly, Q used the trans- 
formation of coordinates, x = X2 cos(0i) + 1/2 sin((/>i) and 
y = —X2 sin(c/)i) + j/2 cos(</>i), and plotted the limit cycle 
in the (x, y) plane. Then the trajectory on this trans- 
formed plane looks like a ring. 

Below the boundary of phase synchronization regime 
(when £ < £ c ), the dynamics of the phase difference 9(t) 
is persistently and intermittently interrupted by sudden 
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phase slips (2ir jumps). As e increases further, the limit 
cycle starts to enclose the origin near the first critical 
point and the phase destruction begins. Two critical val- 
ues e t = 0.1097 and e c = 0.124 were observed @, such 
that (i) the phase synchronization occurs for e > e c , (ii) 
the (ring) intermittent behavior occurs for et < e < e c , 
and (iii) the asynchronous dynamic occurs for e < e t . 
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FIG. 9. (Color online) Ring intermittency in coupled Rossler 
oscillators with parameters ei — 0, £2 = e, and control pa- 
rameters uji — 1.0, L02 = 0.95. The ring intermittency occurs 
between et = 0.1097 and e c = 0.124. (a) Rates n and synchro- 
nization index 7. (b) Relative frequencies of desynchroniza- 
tion events of different durations. Thick solid line stands for 
actual duration, while thin solid line represents the estimates 
from the transition rates. 

Fig.EJa) shows that during the ring intermittency, the 
rate r± slowly decreases. However the rates r^,3 change 
drastically during the ring intermittency while is an 
almost constant during this parameter regime. These 
changes of the rates imply that the probabilities of long 
cycles are getting smaller as the coupling strength in- 
creases (see Fig. EKb)). On the other hand, the probabil- 
ities of short desynchronizations (lasting for one or two 
cycles of oscillations) are getting higher. 

Now let us compare the ring intermittency in Fig. [5] 
with the eyelet intermittency in Fig. [H Both systems 
are unidirectionally coupled Rossler oscillators, but with 
different control parameters Wi,2- 

Since the phase slips in the eyelet intermittency are 
rare, as one can expect, rate n in Fig. [H^a) has much 
smaller value than that in Fig. |9]Ja). Rate r± in Fig. Uta) 
has slightly lower value than that in Fig.Uta). Rate T2 in 
Fig. [Ha) shows huge fluctuation while that in Fig. Uta) 
does not. The probability of long desynchronization 
events (Cycle> 5) for the ring intermittency is close to 
zero while that in the eyelet intermittency starts to de- 
crease from high level (> 0.5). Thus, long desynchro- 
nization events are dominant for the eyelet intermittency 



while short desynchronization events are dominant for 
the ring intermittency. 

V. COUPLED NEURONS EXAMPLE 

Intermittency and intermittent and otherwise imper- 
fect synchronization can be observed in the brain in 
different conditions (for example, see [l?], EH [H, EH). 
Thus it will be useful to apply the method to more re- 
alistic neural-like model systems to study what rates rj 
may reveal about synchronized dynamics. We consider 
two mutually coupled neurons described by a single- 
compartment conductance-based Hodgkin-Huxlcy type 
equations. We couple these neurons through inhibitory 
synapses. The details for this kind of a model can be 
found in, for example, (24[. The equations for each cell 
can be written as: 

C m v = -I L - I Na - I K - I syn + W + J , (13) 
x' = (p(a x (v)(l - x) - /3 x (v)x), 

where I L = g L (v - E L ), I Na = g Na m 3 (v)h(v - E Na ), 
Ik = gKn 4 (v — Ek) represent leak, sodium and potas- 
sium currents respectively. Independent Gaussian white 
noise W with zero mean and standard deviation of 2.12 
is added to both neurons. Iq is an external current, x 
stands for three different gating variables m, n and h with 
the following a(v) and j3(v): 

a n (v) = 0.01(v + 55)/(l - exp(-(v + 55)/10)), (14) 
P n {v) = 0.125 exp(-(v + 65)/80), 
a m (v) = 0.1(t> + 40)/(l - cxp(-(v + 40)/10)), 
p m (v) =4exp(-(v + 65)/18), 
a h (v) = 0.07cxp(-(« + 65)/20), 
p h (v) = l/(l + exp(-(« + 35)/10)). 

The term I syn represents the synaptic current. For a 
cell i e {1,2}, the synaptic current I syn ,i = 9syn,i{vi — 
v syn )sj where j € {1,2} \ {i}. The synaptic variable s 
(the fraction of activated synaptic channels) is modeled 
by the first-order kinetic equation in the form: 

s' = a(l- s)H oo {v-0 v ) -/3s, (15) 

where H^ix) = 1/(1 + exp[— x/a s ]) is a sigmoidal func- 
tion. The parameter values are the following: C' m = 1, 
g Na = 120, E Na = 50, g K = 36, E K = -77, g K = 0.3, 
E L = -54.4, v syn = -85, Io = 10, tp = 0.35, 9 V = 0, 
a s = 5, a = 2, and (3 = 0.05. Both cells are self- 
oscillators in the absence of coupling {g syn .i = 0). Volt- 
ages v for both neurons exhibit spiky time-series (se- 
quence of action potentials) and were filtered to the 
gamma-band (30 — 80 Hz). The Hilbert transforma- 
tion was used to define the phases and then first-return 
maps were constructed. Fig. [10] shows that the synchro- 
nization index 7 changes slowly for a range of values of 
gsyn,2 (for fixed g syn ,i = 0.1), but the rates r» exhibit 
more substantial variations depending on g syn ,2- As can 
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FIG. 10. (Color online) Rates n and synchronization index 7 
for varied g sy n,2 with fixed g syn ,i = 0.1. Although the levels 
of the synchronization index 7 are similar, the rate r2,4 expe- 
rience more substantial variation depending on the synaptic 
strength g ayn ,2- Small noise also exerts some influence on the 
rates, which is not as prominent in 7. Thick solid line is the 
noiseless case and thin solid line is for the noise W . 

be seen from Fig. ITW a. d), the rates r2,4 are changing 
substantially with the similar 7. This implies that the 
rates describe the fine temporal structure of synchroniza- 
tion/desynchronization events. The rates (and thus the 
fine temporal structure of synchronous/desynchronous 
dynamics) may be more affected by noise than synchro- 
nization index 7 (Fig. fTU|) . Especially for small values of 
coupling g sy n,2 overall degree of phase locking is not influ- 
enced by small noise. However the rates V2,3A experience 
substantial variations and thus corresponding timing of 
synchronization/desynchronization events is dependent 
on the noise level. Therefore, exploring the fine tempo- 
ral structures of coupled oscillators in experimental data 
may provide some help in constraining multiple alterna- 
tives in modeling studies, like in [l9j |. 



VI. DISCUSSION 

We studied the first-return maps of the phase- 
difference between two oscillatory signals. If there is 
some phase locking is present on the average, we can 
study the phase-locking relationship on each cycle of os- 
cillations. The phase space of this map is partitioned into 
four regions and the transition rates r, between regions 
are computed. These transition rates characterize the de- 
tails of temporal dynamics, which arc missed by standard 
synchronization measures (such as 7 in Eq. |[T])). 

Using a series of characteristic model systems, such as 
coupled skew tent maps (for which the values of Lya- 



punov exponents can be computed analytically), coupled 
Rossler and Lorenz oscillators, and coupled Hodgkin- 
Huxley neuronal models, we showed that the aforemen- 
tioned transition rates provide a complimentary descrip- 
tion of synchronized dynamics. They describe the dy- 
namics away from the synchronized state, so that even 
if the stability (transversal Lyapunov exponent) of syn- 
chronization manifold and time-series based measure of 
synchrony (such as a phase-locking index) are the same, 
the rates will capture the details of desynchronization 
events. In an extreme case, the same level of synchrony 
may be reached with a relatively rare but long desyn- 
chronization events and numerous short desynchroniza- 
tion events. The rates allow for discrimination of this 
alternative. 

The duration of desynchronization events can be com- 
puted from this approach and is related to the rates 
(although not necessarily defined by them completely). 
The traditional analysis of intermittent behavior, in- 
cluding intcrmittency near synchronization onset, con- 
siders duration of laminar (i.e. synchronized in the 
case of synchronization) episodes. These distributions 
are universal and depend on the type of intermittency. 
But the same intermittency type in different systems 
can naturally have different dynamics of synchroniza- 
tion/desynchronization events (as reinjection mechanism 
is not defined by the type of intcrmittency). We 
showed how the rates are able to capture this difference 
and corresponding difference in the timing of synchro- 
nized/desynchronized dynamics. 

/From the phase-spacc-based approach point of view, 
the system spends a substantial fraction of time away 
from the synchronization manifold, in other parts of the 
phase space. Therefore the stability of this manifold and 
the bifurcations that it can experience provide incom- 
plete information about the dynamics. The introduced 
rates n can help to fill this gap. In the time-series based 
approach, the global synchronization measures are not 
designed to describe the relationship between phases on 
each cycle of oscillations. The approach considered here 
allows to inspect the phase locking at each cycle of os- 
cillations (naturally, only if some synchronization level is 
present overall). Recent experimental and modeling re- 
sults [HI, [l!| indicate that the fine temporal structure 
of synchronization is important. Here we used standard 
test systems to explain how the analysis of this structure 
works. 

The methodology considered here can be applied to 
any appropriate data with some level of phase locking 
regardless of the mechanisms of synchronous dynamics. 
We clarified here what this method yields in terms of 
characterization of the phase space and in terms of the 
fine temporal structure of phase locking. However, we 
think in future studies the method maybe used as a tool 
to study some mechanisms of synchrony or at least distin- 
guish between different mechanisms of synchrony (in par- 
ticular when parameters of coupling cannot be changed 
in experiment). 
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We would like to conclude the paper with four notes re- 
garding the method. First, we want to reiterate that for 
this kind of analysis the signals should be phase-locked on 
average. The synchrony is essentially non-instantaneous 
phenomenon. Only if there is a preferred phase-locking 
angle, we can follow deviations from it on each cycle (and 
that is yet another reason for necessity of computations of 
the global synchronization measures or synchrony mea- 
sure over moderate size time- window, like in (l6| ) . 

Second, the partitioning of the phase space of the <f> x ^ 
map into four parts is somewhat arbitrary. It simplifies 
the computation of rates, as there are only four of them. 
It also implies that if the phases do not deviate from 
more than 7r/2, they are considered in synchrony. This 
appears to be a reasonable assumption for some systems 
(for example, smaller tolerance maybe easily destroyed by 
noise and fluctuations, larger tolerance is less likely to be 
acceptable for information transmission in the synchro- 
nized system). However in certain applications a finer 
partitioning may be necessary. This is apparently pos- 
sible to implement, although the rates will be harder to 
interpret. Moreover, the major idea of the considered ap- 
proach - to study how the phase difference deviates from 
the preferred phase-locking angle - can be considered not 
only in a discrete-time framework as we did here, but 
also can be generalized to continuous time. 



Third, the approach developed here is agnostic to the 
presence or absence of noise in the time-series. Noise 
surely can contribute to the generation of imperfect phase 
locking and affect the transition rates (as in the example 
in the previous Section). The properties of noise may 
be reflected in the transition rates and the considered 
methodology may be a tool to study the effect of noise on 
synchronous dynamics. This appears to be an important 
subject for future research. 

Fourth, if the system is very close to synchronization 
threshold, this analysis may not necessarily be very infor- 
mative, because the system will spend most of the time 
in the synchronized state. However, as we mentioned in 
the Introduction, synchronization in living systems may 
be highly variable. In this case, synchronized events may 
be relatively rare and/or relatively short (although they 
may still be significant functionally). Thus the approach 
considered here (as well as its variations) have a poten- 
tial to describe functionally important temporal features 
of synchronization/desynchronization dynamics. 
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